A Python application for visualizing an imbricate thrust system: Palomeque Duplex (SE, Spain)¶

Al final no hemos hecho esto

We import the math libraries, data analysis libraries and image libraries.

In [1]:
# math libraries

import numpy as np
import sympy

from sympy import *
t = Symbol('t')

from scipy.interpolate import interp2d
from scipy.interpolate import griddata

import math

# data analysis library


import pandas as pd


# import plotting libraries

import plotly.offline as go_offline
import plotly.graph_objects as go
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.ticker as ticker

# imaging library

from PIL import Image

We set the data and imagen directories

In [2]:
DATADIR='data/' # Directory with the data
FIGURESDIR='figures/' # Figures produced

Topography¶

To obtain the topography of the study area we delimited the area and then used the Copernicus DEM tiff image eu_dem_v11_E30N10.Tiff. Since this image weighs 1.1 G, we crop it to get a smaller image that looks like

In [3]:
Im=Image.open(FIGURESDIR+'eu_dem_v11_E30N10_crop.png')
In [4]:
Im
Out[4]:

We build a mesh, each point of this mesh is given by its UTM_X, UTM_Y and elevation. The data for this mesh is saved as a cvs file that we open.

In [5]:
# Topografy data

topo=pd.read_csv(DATADIR+'topografy.csv')
In [6]:
topo
Out[6]:
UTM_X UTM_Y elevation
0 618131.414840 4.197675e+06 702.294556
1 618132.694676 4.197586e+06 707.380493
2 618133.974489 4.197497e+06 722.496460
3 618135.254279 4.197409e+06 721.459595
4 618136.534045 4.197320e+06 732.357910
... ... ... ...
2585 623304.834524 4.195087e+06 1092.907593
2586 623306.169278 4.194999e+06 1152.172607
2587 623307.504009 4.194910e+06 1132.570068
2588 623308.838716 4.194821e+06 1062.827759
2589 623310.173398 4.194732e+06 986.483154

2590 rows × 3 columns

We split the data in topo to get three numpy arrays

In [7]:
topo1=topo[['UTM_X','UTM_Y','elevation']].to_numpy()

topox=[x[0] for x in topo1]
topoy=[x[1] for x in topo1]
topoelv=[x[2] for x in topo1]

In the list cotasxy we save the dimensions of the model.

In [8]:
cotasxy=[min(topox)+200,max(topox)-2000,min(topoy)+100,max(topoy)-200]

The following function surface interpolates the data on a grid of points, the interpolation method can be linear, cubic or nearest.

In [9]:
def surface(grid,metodo): #metodo='linear', 'cubic' o 'nearest'

    tx=grid[0]
    ty=grid[1]
    tz=grid[2]

    X = np.linspace(min(tx), max(tx))
    Y = np.linspace(min(ty), max(ty))
    X, Y = np.meshgrid(X, Y)  # 2D grid for interpolation
    Z = griddata((tx,ty),tz,(X,Y), 
                 method=metodo)
    
    return [tx,ty,tz,X,Y,Z]

We choose the linear interpolation method and apply it to the topographic data.

In [10]:
topo_gr=[topox,topoy,topoelv]
topo_linear=surface(topo_gr,'linear')

We are ready to create the figure that realizes the 3D model. We use two color scales 'YlGn' and 'gray' to get more depth in the model.

In [11]:
fig=go.Figure()

# topografía

fig.add_trace(go.Surface(  x=topo_linear[3],y=topo_linear[4],z=topo_linear[5],
                           colorscale='YlGn',
                           opacity = 0.8,
                           name='Topography',
                           legendgroup='topo',
                           showlegend=True,
                           showscale=False))
fig.add_trace(go.Surface(  x=topo_linear[3],y=topo_linear[4],z=topo_linear[5],
                           colorscale='gray',
                           opacity = 0.3,
                           name='Topo',
                           legendgroup='topo',
                           showlegend=False,
                           showscale=False))

camera = dict(up=dict(x=0, y=0, z=1),
              center=dict(x=0, y=0, z=0),
               eye=dict(x=0, y=-1, z=2) )

fig.update_layout( title="Topography Model",
                  #paper_bgcolor = 'black',
                 scene = dict(
                     xaxis=dict(title='UTM_X', 
                                tickfont = dict(size = 10,color = 'black'),
                                title_font_size=10,
                                titlefont_color='black', 
                                range=[cotasxy[0],cotasxy[1]],
                                backgroundcolor='white',
                                color='black',
                                gridcolor='gray'),
                     yaxis=dict(title='UTM_Y',
                                tickfont = dict(size = 10,color = 'black'),
                                title_font_size=10,
                                titlefont_color='black', 
                                range=[cotasxy[2],cotasxy[3]],  
                                backgroundcolor='white',
                                color='black',
                                gridcolor='gray'),
                      zaxis=dict(nticks=4,
                                tickfont = dict(size = 10,color = 'black'),
                                title='Elevation', 
                                title_font_size=10,
                                titlefont_color='black', 
                                range=[500,1000],
                                backgroundcolor='white',
                                color='black', 
                                gridcolor='gray'),
                     aspectratio=dict(x=2, y=1.8, z=0.3)),
                     #aspectmode='data'),
                     scene_camera= camera
                 ) 
              

fig.show()

Contacts and faults¶

With the following function contact_tr we process the data from the cvs files in which the stratigraphic and mechanical data are stored

In [12]:
def contact_tr(csv,aumento_elevacion):
    tr=pd.read_csv(csv)
    tr1=tr[['x','y','SAMPLE_1']].to_numpy()
    tx=[x[0] for x in tr1]
    ty=[x[1] for x in tr1]
    telv=[x[2]+aumento_elevacion for x in tr1]
    return [tx,ty,telv]

Load and process the data

In [13]:
# mechanic contacts ????

palomeque=contact_tr(DATADIR+'palomeque.csv',0)+['palomeque_thrust']
palomeque_ds=contact_tr(DATADIR+'palomeque_ds.csv',0)+['palomeque_desplazado']
calvillo=contact_tr(DATADIR+'calvillo.csv',0)+['calvillo_thrust']
calvillo_ds=contact_tr(DATADIR+'calvillo_ds.csv',0)+['calvillo_desplazado']

# faults

fault2=contact_tr(DATADIR+'fault2.csv',0)+['fault2']
fault1=contact_tr(DATADIR+'fault1.csv',0)+['fault1']
fault3=contact_tr(DATADIR+'fault3.csv',0)+['fault1']

# stratigraphic contacts ?????

E12cal=contact_tr(DATADIR+'E1E2 calvillo.csv',0)+['E1-E2 calvillo']
E12pal=contact_tr(DATADIR+'E1E2 palomeque.csv',0)+['E1-E2 palomeque']
E23pal=contact_tr(DATADIR+'E2E3 palomeque.csv',0)+['E2-E3 palomeque']
E3Opal=contact_tr(DATADIR+'E3O palomeque.csv',0)+['E3-O palomeque']
PEpal=contact_tr(DATADIR+'PE palomeque.csv',0)+['P-E palomeque']
PEcal=contact_tr(DATADIR+'PE calvillo.csv',0)+['P-E calvillo']
E12cal_ds=contact_tr(DATADIR+'E1E2 calvillo_ds.csv',0)+['E1-E2 calvillo ds']
E12pal_ds=contact_tr(DATADIR+'E1E2 palomeque_ds.csv',0)+['E1-E2 palomeque ds']

The following function contact_dat uses the comad Scatter3d to create de contact data for the figure.

In [14]:
def contact_dat(x,y,elv,h,mode,nam,color,gr,t,ls):
    z=[x+h for x in elv]
    return go.Scatter3d(x=x, y=y, z=z,
            mode =mode,
            name=nam,
            legendgroup=gr,
            showlegend=t,
            line=dict(color=color,
                      dash=ls,
                      width=5)
                        )
In [15]:
cont=[palomeque_ds,palomeque,calvillo,calvillo_ds]
cab=[E12cal,E12pal,E23pal,E3Opal,PEpal,PEcal,E12cal_ds,E12pal_ds]
fal=[fault1,fault2,fault3]

contacts=[contact_dat(palomeque_ds[0],palomeque_ds[1],palomeque_ds[2],30,
                        "lines",'Thrusts','black','contacts',True,None)
          ]+[contact_dat(x[0],x[1],x[2],30,"lines",x[3],'black','contacts',False,None) for x in cont[1:]]

cabs=[contact_dat(E12cal[0],E12cal[1],E12cal[2],30,"lines",'Stratigraphic contacts','black',
                            'cabalgamientos',True,'dot')
              ]+[contact_dat(x[0],x[1],x[2],30,"lines",x[3],'black','cabalgamientos',False,'dot') for x in cab[1:]]

faults=[contact_dat(fault1[0],fault1[1],fault1[2],50,"lines",'Strike-slip faults','black','faults',True,None)
      ]+[contact_dat(x[0],x[1],x[2],30,"lines",x[3],'black','faults',False,None) for x in fal[1:]]
In [16]:
fig=go.Figure(contacts+cabs+faults)

# topography

fig.add_trace(go.Surface(  x=topo_linear[3],y=topo_linear[4],z=topo_linear[5],
                           colorscale='YlGn',
                           opacity = 0.8,
                           name='Topography',
                           legendgroup='topo',
                           showlegend=True,
                           showscale=False))
fig.add_trace(go.Surface(  x=topo_linear[3],y=topo_linear[4],z=topo_linear[5],
                           colorscale='gray',
                           opacity = 0.3,
                           name='Topo',
                           legendgroup='topo',
                           showlegend=False,
                           showscale=False))

camera = dict(up=dict(x=0, y=0, z=1),
              center=dict(x=0, y=0, z=0),
               eye=dict(x=0, y=-1, z=2) )

fig.update_layout( title="Topography Model",
                  #paper_bgcolor = 'black',
                 scene = dict(
                     xaxis=dict(title='UTM_X', 
                                tickfont = dict(size = 10,color = 'black'),
                                title_font_size=10,
                                titlefont_color='black', 
                                range=[cotasxy[0],cotasxy[1]],
                                backgroundcolor='white',
                                color='black',
                                gridcolor='gray'),
                     yaxis=dict(title='UTM_Y',
                                tickfont = dict(size = 10,color = 'black'),
                                title_font_size=10,
                                titlefont_color='black', 
                                range=[cotasxy[2],cotasxy[3]],  
                                backgroundcolor='white',
                                color='black',
                                gridcolor='gray'),
                      zaxis=dict(nticks=4,
                                tickfont = dict(size = 10,color = 'black'),
                                title='Elevation', 
                                title_font_size=10,
                                titlefont_color='black', 
                                range=[500,1000],
                                backgroundcolor='white',
                                color='black', 
                                gridcolor='gray'),
                     aspectratio=dict(x=2, y=1.8, z=0.3)),
                     #aspectmode='data'),
                     scene_camera= camera
                 ) 
              

fig.show()

We add texts for heights and Stratigraphic levels to the figure

In [17]:
namesx=[619200,620800]
namesy=[4196500,4196250]
namesz=[850,850]
namestxt=['Calvillo height','Palomeque height','A','A$^\prime$']

tx=[calvillo[0][1]+50,calvillo[0][3],calvillo[0][4],calvillo[0][5],calvillo[0][7]]
ty=[calvillo[1][1]-50,calvillo[1][3]-50,calvillo[1][4]-70,calvillo[1][5]-50,calvillo[1][7]-70]
tz=[calvillo[2][1]+40,calvillo[2][3]+30,calvillo[2][4]+50,calvillo[2][5]+15,calvillo[2][7]+20]
In [18]:
nnx=[618900,619350,620200,619350,619400,620120,620700,
     618850,619200,620900,621100,619700]
nny=[4196800,4196700,4196400,4196200,4196000,4196100,4196000,
     4195500,4195500,4195800,4195500,4195200]
nnz=[700,750,750,750,700,790,700,750,700,600,650,700]
nntxt=['M1','P','P','E1','E2','E1','E2','E1','E2','E3','O2','E1']
In [19]:
fig=go.Figure(contacts+cabs+faults)

# topography

fig.add_trace(go.Surface(  x=topo_linear[3],y=topo_linear[4],z=topo_linear[5],
                           colorscale='YlGn',
                           opacity = 0.8,
                           name='Topography',
                           legendgroup='topo',
                           showlegend=True,
                           showscale=False))
fig.add_trace(go.Surface(  x=topo_linear[3],y=topo_linear[4],z=topo_linear[5],
                           colorscale='gray',
                           opacity = 0.3,
                           name='Topo',
                           legendgroup='topo',
                           showlegend=False,
                           showscale=False))


trace = go.Scatter3d(
                   x=namesx, y=namesy, z=namesz,
                   text=namestxt,
                   name='Heights',
                   mode="text",
                   textfont=dict(color=["black","black"],size=13),
                   hoverinfo="skip")

fig.add_trace(trace)

trace2 = go.Scatter3d(
                   x=nnx, y=nny, z=nnz,
                   text=nntxt,
                   name='Stratigraphic levels',
                   mode="text",
                   legendgroup='cabalgamientos',
                   showlegend=False,
                   textfont=dict(color='black',size=18),
                   hoverinfo="skip")

fig.add_trace(trace2)


camera = dict(up=dict(x=0, y=0, z=1),
              center=dict(x=0, y=0, z=0),
               eye=dict(x=0, y=-1, z=2) )

fig.update_layout( title="Topography Model",
                  #paper_bgcolor = 'black',
                 scene = dict(
                     xaxis=dict(title='UTM_X', 
                                tickfont = dict(size = 10,color = 'black'),
                                title_font_size=10,
                                titlefont_color='black', 
                                range=[cotasxy[0],cotasxy[1]],
                                backgroundcolor='white',
                                color='black',
                                gridcolor='gray'),
                     yaxis=dict(title='UTM_Y',
                                tickfont = dict(size = 10,color = 'black'),
                                title_font_size=10,
                                titlefont_color='black', 
                                range=[cotasxy[2],cotasxy[3]],  
                                backgroundcolor='white',
                                color='black',
                                gridcolor='gray'),
                      zaxis=dict(nticks=4,
                                tickfont = dict(size = 10,color = 'black'),
                                title='Elevation', 
                                title_font_size=10,
                                titlefont_color='black', 
                                range=[500,1000],
                                backgroundcolor='white',
                                color='black', 
                                gridcolor='gray'),
                     aspectratio=dict(x=2, y=1.8, z=0.3)),
                     #aspectmode='data'),
                     scene_camera= camera
                 ) 
              

fig.show()

To add triangles to the contact tracks we need some mathematical functions

In [20]:
def eq(A,B): # A and B are points
    x1, y1 = A
    x2, y2 = B
    # calculate the distance between the two points
    distance = math.sqrt((x2 - x1)**2 + (y2 - y1)**2)

    # calculate the angle between the two points
    angle = math.atan2(y2 - y1, x2 - x1)

    # calculate the coordinates of the third point (the first one)
    x3_1 = x1 + distance * math.cos(angle - (1 * math.pi / 3))
    y3_1 = y1 + distance * math.sin(angle - (1 * math.pi / 3))

    # calculate the coordinates of the third point (the second one)
    x3_2 = x1 + distance * math.cos(angle + (1 * math.pi / 3))
    y3_2 = y1 + distance * math.sin(angle + (1 * math.pi / 3))
    return[[x3_1,y3_1],[x3_2,y3_2]]

    

def tr(i,contacto,h,s,hh):
    A=(contacto[0][i],contacto[1][i])
    c1z=contacto[2][i]+h
    c2x=contacto[0][i+1]
    c2y=contacto[1][i+1]
    d=(c2x-A[0],c2y-A[1])
    nd=np.sqrt(d[0]**2+d[1]**2)
    dn=(d[0]/nd,d[1]/nd)
    B=(A[0]+dn[0]*150,A[1]+dn[1]*150)
    C=eq(A, B)[s]
    return [[A[0],B[0],C[0]],[A[1],B[1],C[1]],[c1z,c1z,c1z+hh]]
In [21]:
tcalvillo=[tr(1,calvillo,15,0,100),tr(3,calvillo,15,0,50),tr(5,calvillo,15,0,50),tr(7,calvillo,15,0,0)]

tcalvillo_dat=[go.Mesh3d(x=x[0],y=x[1],z=x[2],
                        alphahull=5, opacity=1, color='black',
                        i = np.array([0]),
                        j = np.array([1]),
                        k = np.array([2]),
                         legendgroup='contacts',
                         showlegend=False,
                       ) for x in tcalvillo]

tcalvillo_ds=[tr(1,calvillo_ds,25,1,30)]
tcalvillo_ds_dat=[go.Mesh3d(x=x[0],y=x[1],z=x[2],
                        alphahull=5, opacity=1, color='black',
                        i = np.array([0]),
                        j = np.array([1]),
                        k = np.array([2]),
                            legendgroup='contacts'
                       ) for x in tcalvillo_ds]

tpalomeque=[tr(1,palomeque,15,0,70),tr(4,palomeque,15,0,70),tr(5,palomeque,15,0,90),
            tr(6,palomeque,15,0,70),tr(7,palomeque,15,0,30)]

tpalomeque_dat=[go.Mesh3d(x=x[0],y=x[1],z=x[2],
                        alphahull=5, opacity=1, color='black',
                        i = np.array([0]),
                        j = np.array([1]),
                        k = np.array([2]),
                          legendgroup='contacts'
                       ) for x in tpalomeque]

tpalomeque_ds=[tr(1,palomeque_ds,5,1,0)]

tpalomeque_ds_dat=[go.Mesh3d(x=x[0],y=x[1],z=x[2],
                        alphahull=5, opacity=1, color='black',
                        i = np.array([0]),
                        j = np.array([1]),
                        k = np.array([2]), 
                        legendgroup='contacts',
                       ) for x in tpalomeque_ds]

The triangles data are stored in tcont.

In [22]:
tcont=tcalvillo_dat+tcalvillo_ds_dat+tpalomeque_dat+tpalomeque_ds_dat
In [23]:
fig=go.Figure(contacts+cabs+faults+tcont)

# topography

fig.add_trace(go.Surface(  x=topo_linear[3],y=topo_linear[4],z=topo_linear[5],
                           colorscale='YlGn',
                           opacity = 0.8,
                           name='Topography',
                           legendgroup='topo',
                           showlegend=True,
                           showscale=False))
fig.add_trace(go.Surface(  x=topo_linear[3],y=topo_linear[4],z=topo_linear[5],
                           colorscale='gray',
                           opacity = 0.3,
                           name='Topo',
                           legendgroup='topo',
                           showlegend=False,
                           showscale=False))


trace = go.Scatter3d(
                   x=namesx, y=namesy, z=namesz,
                   text=namestxt,
                   name='Heights',
                   mode="text",
                   textfont=dict(color=["black","black"],size=13),
                   hoverinfo="skip")

fig.add_trace(trace)

trace2 = go.Scatter3d(
                   x=nnx, y=nny, z=nnz,
                   text=nntxt,
                   name='Stratigraphic levels',
                   mode="text",
                   legendgroup='cabalgamientos',
                   showlegend=False,
                   textfont=dict(color='black',size=18),
                   hoverinfo="skip")

fig.add_trace(trace2)


camera = dict(up=dict(x=0, y=0, z=1),
              center=dict(x=0, y=0, z=0),
               eye=dict(x=0, y=-1, z=2) )

fig.update_layout( title="Topography Model",
                  #paper_bgcolor = 'black',
                 scene = dict(
                     xaxis=dict(title='UTM_X', 
                                tickfont = dict(size = 10,color = 'black'),
                                title_font_size=10,
                                titlefont_color='black', 
                                range=[cotasxy[0],cotasxy[1]],
                                backgroundcolor='white',
                                color='black',
                                gridcolor='gray'),
                     yaxis=dict(title='UTM_Y',
                                tickfont = dict(size = 10,color = 'black'),
                                title_font_size=10,
                                titlefont_color='black', 
                                range=[cotasxy[2],cotasxy[3]],  
                                backgroundcolor='white',
                                color='black',
                                gridcolor='gray'),
                      zaxis=dict(nticks=4,
                                tickfont = dict(size = 10,color = 'black'),
                                title='Elevation', 
                                title_font_size=10,
                                titlefont_color='black', 
                                range=[500,1000],
                                backgroundcolor='white',
                                color='black', 
                                gridcolor='gray'),
                     aspectratio=dict(x=2, y=1.8, z=0.3)),
                     #aspectmode='data'),
                     scene_camera= camera
                 ) 
              

fig.show()
In [ ]: